% Define approximate molar mass of components
molmass_casamino = 135.5; % g/L

% Total carbon molar mass
molmass_NZCasePlus = 125/molmass_casamino;
clc;
close all;
%% NZ Case Plus 37C
tempmat = monod_mat_NZCasePlus_25C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1);
    
    
hold on;


smat = size(tempmat);


set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'b', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% NZ Case Plus 37C
tempmat = monod_mat_NZCasePlus_30C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

hold on;
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1);
    
    
hold on;


smat = size(tempmat);


set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'k', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% NZ Case Plus 37C
tempmat = monod_mat_NZCasePlus_37C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

hold on;
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1);
    
    
hold on;


smat = size(tempmat);


set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'r', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);
